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ABSTRACT 


(Distribution  Limitation  Statement  No.  2) 


The  response  of  a  layered  half-space  to  a  progressing 
normal  surface  pressure  is  evaluated.  The  half-space  con¬ 
sists  of  a  single  layer  of  a  locking  material  that  acts 
elastically  after  compaction,  and  an  underlying  elastic 
material.  The  surface  pressure  moves  with  a  constant  velocity 
V,  which  is  subseismic  with  respect  to  the  speed  of  wave 
propagation  in  both  the  upper  layer  after  compaction  and  the 
underlying  half-space.  It  is  assumed  that  a  steady-state 
exists  with  respect  to  a  coordinate  axis  attached  to  the 
moving  load. 

The  essentially  subseismic  layer  -  subseismic  half-space 
geometry  leads  to  an  elliptic  problem  that  is  solved  by  a 
finite  difference  iterative  technique.  A  computer  program 
for  evaluating  stresses  and  velocities  at  points  in  the  medium 
is  available  and  results  are  presented  for  a  typical  con¬ 
figuration  of  interest. 
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SECTION  I 


INTRODUCTION 

This  paper  presents  a  study  of  the  problem  of  a 

progressing  normal  pressure  on  the  surface  of  3  layered 

half-space.  The  half-space  consists  of  a  single  layer  of 

a  locking  material  that  acts  elastically  after  compaction 

and  an  underlying  elastic  material.  The  surface  pressure 

moves  with  a  constant  velocity  V,  which  is  subseismic  with 

respect  to  the  speed  of  wave  propagation  in  both  the  upper 

layer  after  compaction  (c  >  c  >  V)  and  the  underlying 

F1  bl 

half-space  (c  >  c  >  V).  It  is  assumed  that  the  load 
2  b  2 

has  been  moving  steadily  for  a  long  time  so  that  a  steady- 
state  (of  plane  strain)  exists  with  respect  to  a  coordinate 
system  attached  to  the  moving  load.  The  theory  is  developed 
for  two  types  of  progressing  pressure  loadings;  (a)  a  case 
in  which  the  progressing  pressure  is  a  step  function  in  the 
coordinate  £  *  x  +  Vt  over  a  range  of  £  very  large,  and  then 
decays  to  zero  at  £  =  00 ,  and  (b)  a  case  in  which  the 
progressing  pressure  has  an  exponential  decay  in  the 
coordinate  £  =  x  +  Vt,  Fig.  1. 

It  should  be  noted  that  the  material  description  for 
the  layer  assumes  a  locking  material  that,  after  compaction, 
acts  elastically.  This  mathematical  model  is  an  approximation 


1 


(a ) 


- 

•  x  o 

/ 

1  ' 

/ — Pe^ 

■ — — X 

z 

— 

L 

(b) 


.  FIG.  I  -  LOADING  CONDITIONS 


FIG.  2- APPROXIMATION  OF  MATERIAL  BY 
LOCKING -ELASTIC  MODEL 
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to  certain  types  of  xsal  materials  as  shown  in  Fig.  2.  The 
locking-elastic  material  has  been  previously  studied  by  Bleich, 
Ref.  '[I]. 


In  a  previous  paper.  Ref.  [2.]:,  the  steady-state  problem 
for  a  superseismic  layer  -  subseismic  half-space  configuration 
was  studied.  For  that  case,  the  problem  was  of  a  mixed'  type, 
(a)  hyperbolic  in  the  superseismic  layer  and  (b;)  elliptic 
in  the  underlying  half-space.  Consequently,  signals  from  the 
underlying  half-space  outran  the  advancing  surface  pressure 
and  the  interaction  stresses  between  the  layers,  i.e.,  d 

zz 

and  0  extended  over  the  entire  plane  -«=  <  £  <  °°.  The 
zi, 

entire  layer  and  half-space  were  stressed,  including  that 
part  of  the  region  which  was  ahead  of  the  moving  load*.  Sharp 
shock  fronts  (F  and  S)  were  present  in  the  superseismic  layer, 
only  behind  the  leading  front  of  the  moving  load,  Fig.  3. 


The  present  problem  differs  from  the  one  treated  in 
Ref.  [2j  in  a  significant  way,  which  is  best  understood  by 
considering  the  layer  to  consist  of  a  nonlinear  hardening 
material,  as  indicated  in  Fig.  2.  While  the  underlying  half- 
space  is  subseismic  in  both  the  present  case  and  in  Ref.  [2], 
the  situation  in  the  layer  is  superseismic  at  low  stress 
level:,  but  subseismic  at  sufficiently  high  stress  levels. 

In  the  present  problem  which  considers  the  nonlinear  material 
and  a  high  stress  level,  there  will,  therefore,  be  a  leading 
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shock  corresponding  to  the  first  P-front  in  Fig.  3,  but  there 
will  be  no  further  shocks,  to  the  right  of  (behind)  this 
leading  front.  Ahead  of  this  shock  front,  the  stresses  are 
inherently  low.  In  the  case  of  the  locking  material,  Fig.  2, 
the  shock  corresponding  to  the  first  P-front  in  Fig.  3  is  now 
identified  with  a  compaction  front.  Fig.  4.  The  "low"  stresses 
ahead  of  the  front  are  assumed  to  vanish  to  permit  the  use  of 
the  locking  concept.  It  will  be  further  assumed  that  this 
compaction  front  is  plane.  This  assumption  is  correct 
according  to  Ref.  [1]  in  an  infinite  locking-elastic  medium 
subjected  to  a  progressing  step  wave,  but  qot  for  a  decaying 
pressure  pulse  where  it  is  curved  as  indicated  in  Fig.  4.  The 
curvature  of  this  front  is,  however,  not  important  if  the 
layer  thickness  is  not  too  large  compared  to  the  distance 
describing  the  decay  of  the  applied  surface  pressure.  The 
inclination  of  the  locking  front  in  the  present  analysis  will 
be  selected  to  be  equal  to  the  one  found  in  Ref.  [1]  for  a 
step  pressure  of  intensity  P  =  Pq  .  This  leads  to  the  con¬ 
figuration  shown  in  Fig.  5.  The  problem  becomes  one  of  an 
irregularly  shaped  half-space  with  two  layers,  both  of  which 
are  subseismic  with  respect  to  the  moving  velocity  V,  i.e., 

V  is  smaller  than  the  velocity  of  dilatational  (P)  and 
equivoluminal  (S)  waves,  in  both  layers.  Hence,  the  dif¬ 
ferential  equations  of  motion  in  both  layers  become  elliptic 

There  is  a  transition  level  at  which  the  stress  is  in¬ 
sufficiently  high  and  shear  shock  may  occur.  It  is 
assumed  in  the  present  problem  that  the  stress  level  is 
sufficiently  high  so  that  this  does  not  occur. 
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STRESSED  REGION  C  >  C*  >V 

OUTRUNNING  WAVES  2  2 


FIG.  3  -  SUPERSEISMIC  LAYER  —  SUBSEISMIC  HALF  SPACE 
GEOMETRY  OF  REFERENCE  [2] 


MATERIAL  Cp  >  Cs  >  V 

STRESSED  HALF  SPACE  2  2 

FIG. 4  -  GEOMETRY  OF  PRESENT  PROBLEM 
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and  the  solution  reduces  to  that  of  an  elliptical  boundary 
value  problem  in  the  irregular  domain  of  Fig.  5. 

The  elliptic  nature  of  the  problem  leads  to  a  system  of 
simultaneous  linear  algebraic  equations  on  the  displacements 
at  each  point  of  a  finite  difference  grid  that  covers  the 
upper  lay.er  and  the  underlying  half-space.  A  crucial  approxi¬ 
mate  procedure,  which  is  used  in  the  numerical  solution  of  the 
boundary  value  problem,  should  be  discussed  at  this  point. 

The  actual  domain  of  the  boundary  value  problem  discussed  above 
extends  to  infinite  in  three  directions,  i.e.,  -®  <  £  <  00  and 
0  £  z  <  °°,  Fig.  6.  To  obtain  a  finite  number  of  algebraic 
simultaneous  equations  in  the  numerical  analysis,  the  infinite 
domain  is  arbitrarily  reduced  to  the  finite  domain  ABCDEF  as 
shown  in  Fig.  6,  by  assuming  that  the  field  quantities  along 
the  boundaries  of  the  rectangle  BFED  are  same  as  those  that 
would  be  produced  if  a  half-space  of  the  underlying  material 
was  subjected  to  the  same  progressing  pressure  signal  as  in 
tlie  actual  problem,  Fig.  7.  Hence,  it  is  assumed  that  the 
values  of  the  field  quantities  along  the  boundary  BCDEF  of 
Fig.  6  will  be  those  which  occur  in  the  same  locations  for 
the  geometry  of  Fig,  7.  These  field  quantities  may  be  obtained 
by  a  direct  integration  of  the  Cole-Huth  solution.  Ref.  [3]. 
in  a  similar  manner,  the  field  quantities  on  the  line  AB  are 
approximated  by  values  from  the  locking-elastic  media  study 
of  Ref.  [ 1]  for  the  step  pressure  applied  on  an  infinite  half- 
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ASSUMED  STRESSLESC  REGION^N^-  PLANE  COMPACTION  FRONT 

^STRESS  FREE  BOUNDARY  Cp^Cs^V  ( SUBSEISMIC  -ELASTIC) 

Cp^>C<;>V  (SUBSEISMIC  —  ELASTIC  HALF  SPACE) 

5>z 

FIG.  5 -APPROXIMATE.  GEOMETRY  FOR  PRESENT  PROBLEM 


FIG.  6-*  BOUNDARY  CONDITIONS  AT  LARGE  £  AND  Z 
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FIG.  8o- PLANE  APPROXIMATION  FOR  COMPACTION  FRONT  IN  LAYER 


space  and  by  an  integration  of  the  Cole-Huth  solution  for  a 
decaying  pressure. 

In  defense  of  the  approximation  at  the  boundaries,  it  is 
noted  that  it  is  inherent  that  the  stresses  at  far  away 
boundaries  do  not  effect  the  solution  in  the  vicinity  of  the 
shock,  i.e.,  the  only  region  in  which  a  steady-state  solution 
has  actual  meaning.  A  computer  . code  for  the  solution  of  the 
problem  in  the  manner  outlined  has  been  developed  and  is 
available.  Numerical  results  are  shown  for  a  typical  con¬ 
figuration  of  interest. 

Section  II  presents  the  equations  for  the  basic  formu¬ 
lation  of  the  problem.  The  method  of  solution  and  some 
comments  on  the  computer  code  are  given  in  Section  III. 
Finally,  Section  IV  presents  some  numerical  results  and  con¬ 
clusions  . 

This  paper  is  one  of  a  series  of  steady-state  solutions 
that  have  been  obtained  for  various  elastic  and  inelastic 
materials,  Refs.  [4]-[8J.  While  the  solutions  have  been  of 
interest  in  themselves,  an  additional  purpose  was  to  provide 
check  results  for  several  of  the  large  finitt  iifference  and 
finite  element  computer  codes,  which  are  presently  being  used 
to  study  ground  shock  problems  in  elastic  and  inelastic  media. 
Although  these  codes  are  primarily  used  for  the  solution  of 
transient,  rather  than  steady-state  problems,  they  generally 
can  also  be  used  to  model  steady-state  situations  of  the  type 
considered  here. 
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SECTION  II 


FORMULATION  OF  THE  PROBLEM 

Consider  the  geometry  shown  in  Fig.  5  in  which  a  layer 
of  thickness  H  of  a  locking-elastic  solid,  overlies  a  half¬ 
space  of  a  second  elastic  material.  The  layer  is  composed 
of  a  locking  material,  which,  upon  reaching  a  certain  compacting 
strain  ,  thereupon  acts  as  a  linear  elastic  homogeneous  and 
isotropic  material.  Fig.  1.  The  surface  z  =  0  is  subjected  to 
a  normal  surface  pressure  which  moves  with  a  constant  velocity 
V,.  The  velocity  V  of  the  moving  load  is  subseismic  with 
respect  to  dilatational  and  equivoluminal  wave  velocities  both 

in  the  compacted  layer,  i.e.,  c  >  c  >  V,  and  in  the  under- 

F1  S1 

lying  half-space,  c  >  c  >  V.  A  system  of  Cartesian 

2  2 

coordinates  is  used  in  which  z  is  normal  to  the  layer  surface 
and  x  is  parallel  to  the  surface  in  the  direction  of  the 
moving  load.  The  layer  0  £  z  <  H  is  designated  as  region  1 
and  the  half-space  z  >  H  is  designated  as  region  2.  The 
corresponding  elastic  constants  are  indicated  by  the  proper 
subscript . 

For  both  the  compacted  elastic  layer  and  the  underlying 
half-space,  the  elastic  stress-strain  relations  and  the 
equations  of  motion  for  this  plane  strain  problem  are  written 
ad 
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Eliminating  the  stresses  between  these  equations,  the  dis¬ 
placement  equations  of  motion  become 
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3z2 


=  P 


32u 

3t2 
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(A  +  2p) 


32v 

3z2 


+  (A  +  y) 


a2 
3  u 

3x3z 


+  M 


~2 
3  v 

-  2 
3x 


=  P 


32v 

3t2 


(7) 


For  the  steady-state  problem  under  consideration  a  Galilean 
transformation  £  =  x  +  Vt  is  used  to  replace  x  and  t  by 
means  of  the  relations 


_3_  _  _3_ 

3x  "  3?  *  3t  3£ 


The  stress-strain  relations  and  the  equations  of  motion 


(8) 


become 


°»  -  <17^)1(1-v)  ft  +  v  ft1 


°2Z  (l~2v^V  3£  +  (1  v)  3z  ^ 


o  -  cflil  ...  iiii 

0XZ  "  Gln 


(A  +  2y  -  py2)  +  (A  +  y)  +  p  =  o 

3£2  H3Z  3z 2 


2  2  2 
(A  +  2y)  — ~  +  (A  +  y)  -—^7  +  (y  -  pv2)  =  0 


3z" 


3C3z 


Hole  that 


2  _  A  +  2y 

CP  p 


y 

P 


Equations  (12)  and  (13)  may  be  written  in  the  form 

.2. 


A  g^2  +  B  3£3z 


v  +  C  ^  -  0 
3z'Z 


~2  ~  2  _  2 

i>  m  +  ®  -y- +  f  m  -  0 

3z2  3552  3£2 


where 


2  2 
A  =  cp  -  V 


D  =  c 


„  2  2 
B  =  CP  ’  CS 


E  =  E 


C  =  c 


F  "  CS  '  v 
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For  the  subseismic  problem  under  consideration,  the 
coefficients  of  Eqs.  (17)  are  all  positive  and  hence,  the 
system  of  equations,  Eqs.  (15)— (16)  are  elliptic  in  nature 
Consequently,  the  specified  boundary  conditions  on  the 
boundaries  of  the  domain  will  uniquely  determine  the 


solution  of  the  equations  of  motion  in  the  interior  of  the 
domain.  The  boundary  conditions  are  developed  below* 

From  Fig.  4,  it  is  seen  that  boundary  conditions  must 
be  specified  on  the  surface  2=0,  on  the  compaction  front  of 
the  locking-elastic  layer,  on  the  surface  z=H  of  the  half¬ 
space  ahead  of  the  compaction  front  and  on  the  layer  half¬ 
space  interface,  z=H,  behind  the  compaction  front*  The 
requirements  on  the  stress  and  displacement  'quantities  at 
£  =  «>  and  z  =  03  will  be  discussed  in  the  section  on  the 
numerical  computation  of  the  displacements  and  stresses  by 
means  of  a  finite  difference  grid  which  is  superimposed  on 
the  medium. 

The  boundary  conditions  on  the  surface  z=0  are  .given 
by  the  relations 


0 

zz 

(£, 

<— > 
o 
ii 

N 

-PU 

«) 

Or 

SZ 

(€, 

■— N 

o 

II 
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0 

while 

on 

the 

surface  z= 

H  ahead 

of 

the 

compaction  front 

(i . e.  , 

c 

negative ) 

0 

zz 

(S, 

z  =  H) 

= 

0 

a 

Sz 

(?, 

z  =  H) 

- 

0 

The  boundary 

conditions 

on 

the 

layer 

interface  z=H  behind 

the  compaction  front  require  the  continuity  of  stresses 
and  displacements: 


(18) 


(19) 


(20) 

(21) 
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1^(5,  z=H )  =  u2  (C  >  z=H) 
vx(5,  z=R)  =*  v2(C,  z=H) 


Z2,l(5, 

z-H) 

"  °«#2«* 

N 

II 

rc 

z  =  H) 

=  °£z>2(C» 

z=H) 

The  problem  of  a  step  load  moving  on  the  surface  of  a 

half-space  of  a  locking  material  which  upon  compaction, 

becomes  elastic  and  subseismic  ,  i.e.,  c  >  c  >V,  has 

P2  S2 

been  solved  in  Ref.  [1],  It  was  showr.  that  a  stress  dis¬ 
continuity  in  the  form  of  a  compaction  front  occurs  in  the 
material  and  moves  with  the  progressing  surface  pressure. 
For  a  step  wave  surface  pressure  loading,  the  compaction 
front  will  be  plane,  and  the  stress  on  the  front  will  be 
normal  to  the  front  with  no  tangential  component.  The 
normal  strcess  at  the  compaction  front  is  given  by  the 
relation 


°N  ~  ”  v  +  ( 1- 2v) 3  cot  6 

where  3  is  the  angle  of  the  front  with  the  ourface  z=0, 
Fig,,  5a.  The  tangential  stress  on  the  compaction  front  is 
zero. 


(23) 


a 


(24) 


The  material  after  locking  must  have  both  its  dila- 
tational  and  equivoluminal  velocities  greater  than  the 
velocity  V  of  the  progressing  load.  This  is  a  neces¬ 
sary  condition  for  the  continuity  of  the  material 
behind  the  compaction  front,  i.e.,  a  condition  that 
prevents  the  locking  front  from  separating  from  the 
material.  1 
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The  angle  of  inclination  B  of  the  compaction  front  is  given 
in  terms  of  the  velocity  V  of  the  moving  surface  pressure 


where 


y 2  _  _P_  (1-1>) _ _ 

v  sin^B  +  (1— 2v)  sin  B  cos  B 


0  <  B  <  B 

—  —  cr 


where  the  critical  angle  8cr  is  defined  by  the  relation 


tan(2B  ) 
cr 


-  -2 (l-2v) 


(.27) 


The  fact  that  no  steady-state  solution  is  found  if  the 

velocity  V  is  less  than  a  critical  value  V-  corres.ponding 

to  B  is  discussed  in  detail  in.  Section  III  of  Ref.  (l'J 
cr 

and  will  not  be  repeated  here.  It  will  be  assumed  that 
the  B-V  relation  for  the  present  problem  is  such  that  a 
steady-state  solution  exists  and  that  the  locking  front 
meets  the  leading  edge  of  the  moving  pressure  loading  at 
the  point  A  in  Fig.  5a. 


As  discussed  in  Section  I  of  this  paper,  the  com¬ 
paction  front  for  the  present  layered  geometry  is  approxi¬ 
mated  by  a  plane  front  as  shown  in  Fig.  5b.  The  boundary 
conditions  on  this  front  are  given  by  the  relations 


(fr  .  (l-v)P _ 

N^,Z;  "  "  V  +  (l-2v) B  cot  B 


aT(?,z)  -  0 
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where  £  and  z  are  such  that 

*  '  -tan_1$  =  -  |  (30) 

and  the  relation  between  j3  and  V  is  given  by  Eq.  (25). 

The  method  of  solution  of  this  boundary  value  problem 
is, discussed  in  Section  III. 
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SECTION  III 
METHOD  OF  SOLUTION  ^ 

The  equations  of  motion,  Eqs.  (15)- (16)  ar^integrated 
numerically  at  the  pivotal  points  of  a  finite  difference  grid- 
work  which  is  superimposed  on  the  layer  -  underlying  half-s^p-a^ce 
geometry.  Fig.  9.  For  the  elliptic  problem  being  solved ,  these 
differential  equations  are  transformed  into  a  set  of  simultaneous 
linear  algebraic  equations  in  the  displacements  u  and  v.  As 
described  in  Section  I,  this  infinite  system  of  equations  is 
truncated  by  applying  approximate  boundary  conditions  at  all 
pivotal  points  which  lie  on  the  boundaries  of  a  rectangle;  these 
boundaries  are  located  at  a  sufficiently  large  distance  from  the 
area  of  prime  interest,  i.e.,  the  neighborhood  of  the  compaction 
front  in  the  layer  and  the  shallow  portions  of  the  underlying 
half-space.  Once  the  displacements  are  obtained,  the  stresses 
can  be  calculated  at  the  corresponding  grid  points  by  means  of 
finite  difference  operators. 

The  method  of  solution  for  the  system  of  simultaneous 
algebraic  equations  was  dictated  by  the  available  core  size  of 
the  CDC  6600  computer  that  was  utilized.  An  M  x  N  gridwork  of 
pivotal  points  leads  to  an  unsymmetrical  set  of  algebraic 
equations  in  2MN  unknowns.  The  band  width  of  the  equations 
are  of  the  order  of  4N  and  the  core  required  for  a  standard  type 
of  inversion  solution  is  of  the  order  of  8MN  .  For  the  large 


—  FINITE  DIFFERENCE  GRID  AND  AXIS  GEOMETRY 
FOR  NUMERICAL  SOLUTION 


number  of  points  required  to  obtain  acceptable  detail  in  the 
present  problem,  the  direct  inversion  approach  became  im¬ 
practicable.  Consequently,  an  alternative  iterative  technique 
was  used  in  which  the  displacements  u  and  v  at  each  point  in 
the  n+1  approximation  were  calculated  in  terms  of  the  corre- 
sponding  values  in  the  displacement  field  at  the  n  approxi¬ 
mation,  as  required  by  the  appropriate  finite  difference 
operator  at  the  pivotal  point  in  question.  In  this  iterative 
procedure,  only  2MN  of  core  was  required  and  hence,  this 
procedure  was  adopted. 

A  brief  summary  of  the  various  finite  difference 
equations  which  were  utilized  in  the  solution  is  presented 
in  this  section. 


The  partial  differential  equations,  Eqs.  (15)-(16)  may 
be  written  at  all  ordinary  interior  pivotal  points  of  the 
grid  in  terms  of  unequally  spaced  finite  difference  operators. 
Using  the  geometry  of  Fig.  10,  the  expressions  for  the 
various  second  derivative  operators  are  given  by 


n  -  JULtili  L J  I 
5  h£  +  hr 

D  =  ji  ,3+1.1  -  UJ-l] 

z  h  +  h. 

a  b 


_ 2 _ 

hrh£ (hr  +  V 


[hr[i-l,j]  -  (hr  +  h£)[i,j]  + 

+  h4[i+l,j]] 


>  (31) 


(32) 
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FI6.II-  DISPLACEMENT  RESOLUTION  ON  COMPACTION  FRONT 
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BL  ■  h-OhViT)  th  li.j-l)  -  (ha  +  hb)[l,J]  +  haU,j+lJJ  (33) 

aba  b 

c5z  *  ThTTirm  -7- 577  [fi+i.J+U  -  U+1.3-1]  -  ti-i,J+-U  + 

50  r  a  b 

+  (34) 


Substituting  Eqs.  (31)  —  (34)  into  Eqs.  (15)-(16)  leads  to  a 
set  of  two  simultaneous  linear  algebraic  equations  on  the 
values  of  the  displacements  u  and  v: 


bJL(hr  +  h^)  Ui-l,j  "  hrh£  ui,j  + 

+  - § -  (v 

2(hr  +  h£)(ha  +  hb)  VVi+l,j+l 

c  c 

+  ha(ha  +  hb)  Ui,j-1  ~  TTh^  ui,j 

F  F 

ha(ha  +  V  i>i~1  ’  hahb  V±’j 

2(ha  +  bb)(hr  +  h£)  lUi+l,j+l 

D  D 

h£(hr  +  V  Vi_1»j  ~  hrh2,  Vi»J 


hr(hr  +  V  Ui+1>3  + 

Vi+l,j-l  "  Vi-l,j+l 

+  hb(haC+  V  Ui'J+1 

F 

h.(h  +  h  )  vi,j+l 

b  a  d 

ui+l,j-l  "  ui-l,j+l 

D 

hr(hr  +  V  Vi+1’J 


+  Vi.ri5  + 

*  0  (35) 


ui-lrj-l^ 


=  0  (36) 


Special  formulas  are  of  course  required  to  express  the 

/ 

various  conditions  at  the  boundaries  of  the  problem.  A  brief 
description  of  the  boundary  conditions  in  finite  difference 


form  follows. 


li  Boundary  conditions  at  the  compaction  front 


At  the  compaction  front  in  the  locking  material  layer, 

the  normal  stress  is  given  by  Eq.  (23)  and  the  tangential 

stress  is  set  equal  to  zero.  These  conditions  are  written 

in  terms  of  the  displacement  components  and  in  the  £,z  system 

of  coordinates.  Writing  the  stress-displacement  relations  in 

the  normal— tangential  coordinate  system.  Fig.  11, 

3u  3u  3u 

N  3xn  3Xj  3xn 

3un  3ut 
“i  '  +  3Xn] 

and  using  the  N,T  -  £,z  transformation  relations. 


u„  *  -u  sin  B  +  v  cos  3 
N 


uT  =  -u  cos  3  -  v  sin  3 


=  -sin  3 

3xN 


cos  3 


3£  c 

W-  -  -cas  6 


-  -sin  8 

3xt 


the  boundary  conditions  become 


(A  +  y)T  -  (Q  cos  23  +  R  sin  23)  =  af 


Q  sin  23  -  R  cos  23  =  0 
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where 


P  =  ur  +  v 

4  z 

Q  =  VL-  -  V 

4  z 

fc  =  u  4-  v_ 

Z  4 


Simplifying  Eqs.  (40)  to  the  form 


(A  +  y)P  -  gin  26  R  ~  a0 


Q  sin  20  -  R  cos  26-0 


and  substituting  the  finite  difference  operators  of  Eqs-  (31) 
for  the  displacement  derivatives  in  P,  Q  and  R,  Eqs.  (42)  can, 
be  solved  for  the  displacements  u  and  v,.  thus  giving  the 
boundary  conditions  for  the  pivotal  points  which  lie  oh  the 
compaction  front. 


2 .  Boundary  conditions  on  the  surface  z=0  and  on  the  interface 
z=H  (interface  ahead  of  the  shock  front) 


The  boundary  conditions  on  the  surface  z=0  are  given  by 
Eqs.  (18)-(19).  Using  the  finite  difference  operators  of 
Eqs.  (31),  these  conditions  are  written  in  terms  of  the 
boundary  displacements  as 


Ui,j  “  Ui,j  +  1  +  K2(Yi+l,j  "  Vi-l,j) 


v,  .  =  K_(u.,1  ,  -  u.  -  ,)  +  v.  +  KCP 

i,j  3  i+1 , j  i-l,j  i.j+1  5 
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where 

' :  h  K  A  h 

-ir  =  .  .  V  TT  — .  “ _  rr  __  _ P 

2  hr  +  »  3  X  +  2u  *  5  X  +  2y 

■For.  the  stress  free  surface  z=H  which  lies  ahead  of  the  com¬ 
paction  front-,  the  quantity  P  is  set  equal  to  zero  in 
Eqs ..  .(43)- (44). 

3.  Continuity  conditions  on  the  interface  behind  the 
compaction  front  (z=H ,  g  >  0) 


The  continuity  conditions  on  this  interface  are  given 
by  Eqs.  (2-2).  The  equations  on  the  normal  and  tangential 
stresses  at  this  interface  are  given  in  finite  difference 
foTfli  by  the  relations 


A2.  r 

^h,  +  h  ^  i+l,j 


,  .  }l  +  2*1  .  *2  +  2U2, 
Ui-l,j]  +  (  hr +  h7  )vi,j 


*2  +  2^2,  .  A  +  2vl, 

(  h,  )Vi,j+l  (  h  )Vi,j-l 


(46) 


V,  -  V 

(lT  +  \)Ui,j  +  (h  +  h  )(vi+l,j  '  vi-l,j] 


i,j+l 


(47) 


Equations  (46)-(47)  are  solved  for  the  displacements  u.  „  and 

1  *J 

v.  .  ,,  thus  giving  the  displacements  on  the  surface  z=H. 

1 » J 
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A  special  case  occurs  for  the  pivotal  point  which  lies 
both  on  the  compaction  front  and  on  the  interface.  At  this 
point,  one-sided  forward  differences  must  he  used  in  t,he 
expression  for  the  operator  and  £qs»  .(46)—  £47)  become:: 


A..  A  +  2p_  A  -+  2y..  A-  An 

h~  Ui,j  +  (  h^  +  h^  “  (hJ,  +  hx  " 


X2  +  *1  +  V+_2^ 

hA  +  hr  ui~i,  j  (  ha  )vi,3-l  (  \  ^i.Ji-l 


/l  v,  Vi 

h  +  h  Ui,j  h  Vi,J 
a  b  r 


T  U  _  .  , 

ha  r,j-l 


V2 

+  ~ —  u,  _  •+ 
hfa  a  ,3+1 


V. 


+  ( 


hf  +  hr 


h^’i+l.j 


v. 


h„  +  h  i— 1,3 

it  r  * 


Equations  (48)- (49)  may  he  solved  foT  the  displacements 

u.  „  and  v,  a  at  the  point  in  question 
i,J  i,J 


4.  Special  case  for  point  which  is  both  on  the  surface  z—0 
and  the  compaction  front.  Fig.  12 


In  order  to  express  the  requirement  that  the  stress  in 

the  z-direction,  a  ,  at  the  point  Q,  (i,j)  in  fig.  12  be 

z  z 

equal  to  -P,  it  is  necessary  to  define  an  operator  for  3/3z 

at  this  point.  This  is  done  by  utilizing  a  fictitious  point 

M  with  coordinates  (i,j+l)  and  displacements  u^  and  v^  . 

The  boundary  conditions  at  the  surface,  Eqs.  (1.8)—  (19)  are 

written  in  terms  of  the  four  unknown  disp lacement s  u.  .  , 

v.  .  ,  u  and  v  : 

x,j  M  M 


i(48) 


(49) 
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SPECIAL  CASE  FOR  POINT  Q 


(50) 
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0  v  3  t  -  '  ~ 


Approximate  boundary  conditions  at  large  jg|  and  z. 
Figures  6  and  7 


In  Section  I,  the  method  used  to  truncate  the  infinite 
s.atof  finite  difference  equations  on  the  displacements  was 
described;,  Essentially,  this  was  done  by  applying  approxi¬ 
mate  boundary  conditions  at  pivotal  points  on  the  boundaries 
of  a  rectangular  domain  BCDF.F  located  at  large  distances 

and  z)  from  the  domains  of  prime  interest  near  the  com¬ 
paction  front.  These  displacement  quantities  on  the  bounda- 
r'fes  BCDEF  of  the  subseismic  half-space  were  obtained  by 
direct  integration  on  the  computer  of  the  corresponding 
quantities  from  the  Cole-Huth  solution  for  a  point  load  P 
moving  with  subseismic  velocity  on  the  surface  of  a  homo¬ 
geneous  isotropic  elastic  medium,  Ref.  [3].  The  Cole-Huth 
displacements  for  the  point  load  are  given  by  Eqs.  (25)  and 
(26)  of  Ref.  [3]  and  are  repeated  here  for  convenience: 


u 


Btk2 

p  -y  (i 


V 


_p_ 

uu 


[K2  £ n  rT 


6l  Kj  In  rj 


(55) 

(5fc> 


where 


2  -  m: 


ki  “ 


(2  -  M2)2  -  4Bt8l 


2Bt 


K  K 


(2  -  M2)2  -  4$tBl 


(57) 

(58) 
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(59) 


Bl  -  (1  -  li])h  ;  3T  =  (1  - 

and 

(i0L) 

?L  =  X  +  i8LZ  =  rL6 

(ieT) 

?T  =  x  +  i  Bij>  z  =  rTe  j  (60) 

°  i  0L  ;  0T  £  it 

The  distributed  loads  of  the  present  problem  were 
resolved  into  a  set  of  equivalent  concentrated  loads  for 
the  determination  of  these  boundary  conditions.  A  com-r 
prehensive  study  was  made  in  order  to  assess  the  sensitivity 
of  the  displacements  at  pivotal  points  on  the  rectangular 
boundary  BCDEF,  on  this  resolution  procedure.  This  in¬ 
volved  a  study  of  the  influence  in  the  half  plane  of  the 
resolution  of  the  distributed  load  into  n  concentrated  loads 
of  magnitude  P/n  spaced  at  a  distance  of  Ax  units  apart. 

The  results  showed  that  the  resolution  of  the  distributed 
loads  into  n  concentrated  loads  does  not  practically  in¬ 
fluence  the  displacements  of  points  located  at  distances  of 
15  to  20  Ax  from  the  central  loading  point,  in  the  horizontal 
direction,  and  20  to  25  Ax  in  the  vertical  direction.  Hence, 
this  concentrated  load  resolution  approach  was  utilized  for 
the  pres  .it  layered  problem. 
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6 .  Iterative  method  for  the  solution  of  the  system  of  1-Lnear 
algebraic  equations  using  the  over-  and  underrelaxat Ion 
technique 


The  over-  and  underrelaxation  iteration  technique  was 

utilized  in  the  solution  of  the  system  of  simultaneous  algebraic 

equations  on  the  displacements.  This  approach  is  described  in 

detail  in  Ref.  [9]  and  will  be  outlined  in  this  section. 

Assuming  that  the  values  of  the  displacement  field  at  the  n-1 
til 

and  n*"  cycles  are  known,  a  set  of  starting  values  for  the  'Vn+D 
cycle',  u'j  v'  can  be  obtained  by  means  of  a  linear  interpolation 
(under-relaxation)  of  these  previous  values 


u’  =  03u(n)  +  (1-w)  u(n“1) 
v*  «  d)v(n)  +  (l-o;)  v^-15 


w  <  1 


(bl) 


or  a  linear  extrapolation  with  u)  >  1  (overrelaxation).  The 

value  of  the  parameter  w  =  1  would  result  in  using  u^n\  v^n\ 

til 

the  values  of  the  displacement  from  the  n  iterative  cycle 
as  the  starting  values  for  the  n+1  cycle.  Using  the  values 
of  u’  and  v*  ,  Eqs.  (61)  as  starting  values,  the  values  of 
y(n+l)  an(j  y(n+l)  are  0jjtainejj  for  the  n+1  cycle. 


From  the  theory  of  matrix  iterative  analysis,  it  may  be 
shown  that  the  relaxation  factor  w  should  be  chosen  according 
to  the  following  requirements: 

(1)  There  exists  a  limiting  optimum  value  for  which 
the  system  converges  most  rapidly.  If  w  >  ,  the 
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system  converges,  but  with  an  oscillating  pattern 
with  respect  to  successive  values  of  u^n^  and 
a  highly  undesirable  convergence  for  a  digital  com¬ 
puter  since  it  can  lead  to  unacceptable  arithmetic 
conditions.  Consequently,  w  should  he~  chosen  to  be 
smaller  than  this  limiting  value  toy  . 

(2)  The  value  of  depends  on  the  grid  size  and  decreases 
when  the  minimum  spacing  between  two  consecutive 
points  is  decreased. 

(3)  For  large  grids  in  the  present  layered  problem, 

was  found  experimentally  to  be  less  than  unity,  thus 
requiring  the  use  of  underrelaxation,  rather  than 
over relaxat ion . 

(4)  The  value  of  03^  varies  as  u^n^  and  v^n^  approach  their 
respective  true  values  as  the  number  of  iterations  n 
becomes  large.  It  is  therefore  advantageous  to  use 

a  variable  u)^n^  in  the  calculations  in  order  to  speed 
the  rate  of  convergence. 

7 .  Finite  difference  stress  -  displacement  relations 

Once  the  displacements  u  and  v  have  been  obtained  by  the 
iteration,  the  stresses  are  computed  using  the  stress-dis¬ 
placement  relat  ons  in  finite  difference  form: 
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SECTION  IV 


NUMERICAL  RESULTS  AND  DISCUSSION 

This  section  presents  a  set  of  numerical  results  for  the 
case  of  a  progressing  surface  load  with  a  step  distribution 
in  time.  Results  are  presented  for  a  typical  configuration 
that  has  the  following  material  and  geometric  constants? 

c  =  5,000  ft/sec 
1 

\>1  =  0.25 

=  2.85  lb  sec2/ft4 

c  =  8,000  ft/sec 
*2 

v2  =  0.25 

P2  ■  3.80  lb  sec2/ft4 
H  =  100  ft 
PQ  =  1  lb/ft 
V  =  2 , 600  f t/sec 

It  is  of  interest  to  consider  the  typical  variation  of  the 
stresses  O  ,  a,.  and  a^r  a t  points  in  both  the  layer  and 

Z  Z  Z  S  S 

the  half-space,  Figs.  13-15  show  the  variation  of  these 
stresses  with  £  at  the  depths  z  =  60  ft  and  z  -  80  ft  in  the 
layer.  It  is  seen  that  in  each  case,  the  stresses  rise 
sharply  behind  the  shock  front,  oscillate  and  then  decay  to 
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0. 


the  following  limits:  0  T  =  -P.  ,  o  r  T  =  0  ano  ore.  T 

zz  ,L  0  zt,  ,L  t,t,  ,L 

The  peak  value  depends  upon  the  velocity  of  the  moving  load 
and  increases  with  increasing  Mach  numbers  Mp  and  Mg  .  These 
stress  amplifications  are  increased  as  the  transonic  case 
Mg  ®  1  is  approached. 


Figures  16-18  show  the  variation  of  the  stresses  o 


z  z 


anc*  versus  £.  at  points  in  the  underlying  half-space  at 

depth  z  *  260  ft  and  z  =  800  ft.  In  this  case,  the  peak  values 
are  located  on  a  line  that  is  ‘nelined;  the  peak  values  for  the 
deeper  points  occur  for  smaller  values  of  £,  i.e.,  they  run 
ahead  of  the  compaction  front  in  the  layer.  Again,  the  peak 
values  increase  as  Mp  and  I5g  increase  and  approach  high  values 
as  Mg  1,  i.e.,  the  sonic  point  is  approached. 
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